      MODULE PCGUMODULE
        INTEGER,SAVE,POINTER  :: ITER1C,NPC,NORD,NITERC,NNZC,NIAC
        INTEGER,SAVE,POINTER  :: NIAPC,NIWC,NLAPC,NUAPC
        INTEGER,SAVE,POINTER  :: NNZAPROC,NIAPROC
        REAL   ,SAVE,POINTER  :: HCLOSEPCGU,RCLOSEPCGU
        INTEGER,SAVE,POINTER  :: IPRPCGU,MUTPCGU
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   :: LHCHPCGU
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   :: LRCHPCGU
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IT1PCGU
        REAL,             SAVE, POINTER, DIMENSION(:)     :: HCHGPCGU
        REAL,             SAVE, POINTER, DIMENSION(:)     :: RCHGPCGU
        INTEGER,          SAVE, POINTER, DIMENSION(:,:,:) :: NODEC
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:)     :: BC
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:)     :: XC
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:)     :: AC
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:)     :: TAPC
        DOUBLE PRECISION, SAVE, POINTER, DIMENSION(:)     :: APC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IAC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: JAC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IUC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IXMAP
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IAPC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: JAPC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IRAPC
        INTEGER,          SAVE, POINTER, DIMENSION(:)     :: IAPCMAP
C         WORKING ARRAYS        
        INTEGER,         SAVE, POINTER, DIMENSION(:)      :: IWC
        DOUBLEPRECISION, SAVE, POINTER, DIMENSION(:)      :: DC
        DOUBLEPRECISION, SAVE, POINTER, DIMENSION(:)      :: PC
        DOUBLEPRECISION, SAVE, POINTER, DIMENSION(:)      :: QC
        DOUBLEPRECISION, SAVE, POINTER, DIMENSION(:)      :: ZC
        DOUBLEPRECISION, SAVE, POINTER, DIMENSION(:)      :: SCL
      TYPE PCGUTYPE
        INTEGER,POINTER  :: ITER1C,NPC,NORD,NITERC,NNZC,NIAC
        INTEGER,POINTER  :: NIAPC,NIWC,NLAPC,NUAPC
        INTEGER,POINTER  :: NNZAPROC,NIAPROC
        REAL   ,POINTER  :: HCLOSEPCGU,RCLOSEPCGU
        INTEGER,POINTER  :: IPRPCGU,MUTPCGU
        INTEGER,          POINTER, DIMENSION(:,:)   :: LHCHPCGU
        INTEGER,          POINTER, DIMENSION(:,:)   :: LRCHPCGU
        INTEGER,          POINTER, DIMENSION(:)     :: IT1PCGU
        REAL,             POINTER, DIMENSION(:)     :: HCHGPCGU
        REAL,             POINTER, DIMENSION(:)     :: RCHGPCGU
        INTEGER,          POINTER, DIMENSION(:,:,:) :: NODEC
        DOUBLE PRECISION, POINTER, DIMENSION(:)     :: BC
        DOUBLE PRECISION, POINTER, DIMENSION(:)     :: XC
        DOUBLE PRECISION, POINTER, DIMENSION(:)     :: AC
        DOUBLE PRECISION, POINTER, DIMENSION(:)     :: TAPC
        DOUBLE PRECISION, POINTER, DIMENSION(:)     :: APC
        INTEGER,          POINTER, DIMENSION(:)     :: IAC
        INTEGER,          POINTER, DIMENSION(:)     :: JAC
        INTEGER,          POINTER, DIMENSION(:)     :: IUC
        INTEGER,          POINTER, DIMENSION(:)     :: IXMAP
        INTEGER,          POINTER, DIMENSION(:)     :: IAPC
        INTEGER,          POINTER, DIMENSION(:)     :: JAPC
        INTEGER,          POINTER, DIMENSION(:)     :: IRAPC
        INTEGER,          POINTER, DIMENSION(:)     :: IAPCMAP
C         WORKING ARRAYS        
        INTEGER,         POINTER, DIMENSION(:)      :: IWC
        DOUBLEPRECISION, POINTER, DIMENSION(:)      :: DC
        DOUBLEPRECISION, POINTER, DIMENSION(:)      :: PC
        DOUBLEPRECISION, POINTER, DIMENSION(:)      :: QC
        DOUBLEPRECISION, POINTER, DIMENSION(:)      :: ZC
        DOUBLEPRECISION, POINTER, DIMENSION(:)      :: SCL
      END TYPE
      TYPE(PCGUTYPE), SAVE ::PCGUDAT(10)
      END MODULE PCGUMODULE


      SUBROUTINE PCGU7AR(IN,MXITER,IGRID)
C     ******************************************************************
C     ALLOCATE STORAGE FOR PCGU ARRAYS AND READ PCGU DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,   ONLY:IOUT,NCOL,NROW,NLAY,IBOUND
      USE PCGUMODULE
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN)    :: IN
      INTEGER, INTENT(INOUT) :: MXITER
      INTEGER, INTENT(IN)    :: IGRID
C     + + + LOCAL DEFINITIONS + + +
      CHARACTER (LEN=200) :: line
      INTEGER :: i, j, k, n
      INTEGER :: lloc,istart,istop
      INTEGER :: itmem
      REAL    :: r
C     ------------------------------------------------------------------
      ALLOCATE( ITER1C, NPC, NORD, NITERC )
      ALLOCATE( NNZC, NIAC, NIAPC, NIWC )
      ALLOCATE( NLAPC, NUAPC )
      ALLOCATE( NNZAPROC, NIAPROC )
      ALLOCATE( HCLOSEPCGU, RCLOSEPCGU )
      ALLOCATE( IPRPCGU, MUTPCGU )

      ALLOCATE(NODEC(NCOL,NROW,NLAY))
      DO k = 1, NLAY
        DO i = 1, NROW
          DO j = 1, NCOL
            NODEC(j,i,k) = 0
          END DO
        END DO
      END DO
C
C-------PRINT A MESSAGE IDENTIFYING PCG PACKAGE
      WRITE (IOUT,500)
  500 FORMAT (1X,/1X,'PCGU -- UNSTRUCTURED CONJUGATE-GRADIENT SOLUTION',
     &        ' PACKAGE, VERSION 7.01, 08/01/2011',
     &        /1X,8X,'INCLUDES GPU - CUDA AND CPU-MPI SUPPORT')
C
C-------READ AND PRINT COMMENTS, MXITER,ITER1 AND NPCOND
      CALL URDCOM(IN,IOUT,line)
      lloc=1
      CALL URWORD(line,lloc,istart,istop,2,MXITER,r,IOUT,IN)
      CALL URWORD(line,lloc,istart,istop,2,ITER1C,r,IOUT,IN)
      CALL URWORD(line,lloc,istart,istop,2,NPC,r,IOUT,IN)
      CALL URWORD(line,lloc,istart,istop,2,NORD,r,IOUT,IN)
      IF ( NPC.NE.3 .AND. NPC.NE.4 ) NORD = 0
      WRITE (IOUT,510) MXITER, ITER1C, NPC, (NORD>0)
  510 FORMAT (' MAXIMUM OF ',I6,' CALLS OF SOLUTION ROUTINE',/,
     &        ' MAXIMUM OF ',I6,
     &        ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE',/,
     &        ' MATRIX PRECONDITIONING TYPE :',I5,/,
     &        '   NONE:    NPC = 1',/,
     &        '   JACOBI:  NPC = 2',/,
     &        '   ILU0:    NPC = 3',/,
     &        '   MILU0:   NPC = 4',/,
     &        'NEU. POLY.: NPC = 5',/,
     &        ' ILU0 AND MILU0 REORDERING:            ',L5,/,
     &        '   CRS ORDERING:                    NORD = 0',/,
     &        '   SMITH AND ZHANG (2010) ORDERING: NORD > 0')
C
C-------READ HCLOSEPCGU,RCLOSEPCGU,IPRPCGU,MUTPCGU
      READ (IN,*) HCLOSEPCGU,RCLOSEPCGU,IPRPCGU,MUTPCGU
C
C-------PRINT MXITER,ITER1,NPCOND,HCLOSEPCG,RCLOSEPCG,RELAXPCG,NBPOL,IPRPCG,
C-------MUTPCG,DAMPPCG
        WRITE (IOUT,511)
  511   FORMAT (1X,///,36X,'SOLUTION BY THE CONJUGATE-GRADIENT METHOD',
     &        /,35X,43('-'))
      WRITE (IOUT,512) MXITER
  512 FORMAT (1X,19X,'MAXIMUM NUMBER OF CALLS TO PCG ROUTINE =',I9)
      WRITE (IOUT,515) ITER1C
  515 FORMAT (1X,23X,'MAXIMUM ITERATIONS PER CALL TO PCG =',I9)
      WRITE (IOUT,520) NPC
  520 FORMAT (1X,30X,'MATRIX PRECONDITIONING TYPE =',I9)
      WRITE (IOUT,535) HCLOSEPCGU
  535 FORMAT (1X,24X,'HEAD CHANGE CRITERION FOR CLOSURE =',E15.5)
      WRITE (IOUT,540) RCLOSEPCGU
  540 FORMAT (1X,20X,'RESIDUAL CHANGE CRITERION FOR CLOSURE =',E15.5)
      IF ( IPRPCGU.LE.0 ) IPRPCGU = 999
      WRITE (IOUT,545) IPRPCGU, MUTPCGU
  545 FORMAT (1X,10X,
     &        'PCGU HEAD AND RESIDUAL CHANGE PRINTOUT INTERVAL =',
     &        I9,/,1X,4X,
     &        'PRINTING FROM SOLVER IS LIMITED(1) OR SUPPRESSED (>1) =',
     &        I9)
      NITERC = 0
C
C-------ALLOCATE SPACE FOR SOLVER CONVERGENCE ARRAYS
      itmem = MXITER * ITER1C
      ALLOCATE (HCHGPCGU(itmem))
      ALLOCATE (LHCHPCGU(3,itmem))
      ALLOCATE (RCHGPCGU(itmem))
      ALLOCATE (LRCHPCGU(3,itmem))
      ALLOCATE (IT1PCGU(itmem))
C
C-------CALCULATE NUMBER OF ENTRIES IN MODEL GRID
      NNZC  = 0
      NIAC  = 0
      NIAPC = 0
      NIWC  = 0
!      IC    = 0
!      ieq   = 0
      CALL PCGU7CCN(NLAY,NROW,NCOL,NIAC,NNZC,NODEC)
C
C-------ALLOCATE AND INITIALIZE COMPRESSED ROW STORAGE VECTORS
      ALLOCATE(AC(NNZC))
      NIAPC = NNZC
      IF ( NPC.EQ.1 ) THEN
        NIAPC = 1
      ELSE IF ( NPC.EQ.2 ) THEN
        NIAPC = NIAC
      END IF
      ALLOCATE(TAPC(NIAPC))
      ALLOCATE(APC(NIAPC))
      ALLOCATE(IAC(NIAC+1),JAC(NNZC),IUC(NIAC),IXMAP(NIAC))
C       ALLOCATE SPACE FOR ILU0 AND MILU0 NON-ZERO ROW ENTRY VECTOR
      IF ( NPC.EQ.3 .OR. NPC.EQ.4 .OR. NPC.EQ.5 ) THEN
        NIWC = NIAC
      ELSE
        NIWC = 1
      END IF
      ALLOCATE(IWC(NIWC))
C       ALLOCATE SPACE FOR REODERING ILU0 AND MILU0 NON-ZERO ELEMENTS
      NIAPROC  = NIAC
      NNZAPROC = NNZC
      IF ( NORD.LT.1 ) THEN
        NIAPROC  = 1
        NNZAPROC = 1
      END IF
      ALLOCATE(IAPC(NIAPROC+1),JAPC(NNZAPROC),IRAPC(NNZAPROC))
      ALLOCATE(IAPCMAP(NNZAPROC))
C       ALLOCATE WORKING VECTORS FOR PCGU SOLVER      
      ALLOCATE(BC(NIAC),XC(NIAC))
      ALLOCATE(DC(NIAC),PC(NIAC))
      ALLOCATE(QC(NIAC),ZC(NIAC))
      ALLOCATE(SCL(NIAC))
C       INITIALIZE PCGU SOLVER VECTORS
      DO n = 1, NNZC
        AC(n)  = 0.0D0
        JAC(n) = 0
      END DO
      DO n = 1, NIAPC
        APC(n) = 0.0D0
      END DO
      DO n = 1, NIAC+1
        IAC(n) = 0
      END DO
      DO n = 1, NIAC
        IUC(n)   = 0
        IXMAP(n) = 0
        BC(n)    = 0.0D0
        XC(n)    = 0.0D0
C         WORKING ARRAYS
        DC(n)    = 0.0D0
        PC(n)    = 0.0D0
        QC(n)    = 0.0D0
        ZC(n)    = 0.0D0
        SCL(n)   = 0.0D0
      END DO
      DO n = 1, NIWC
        IWC(n)   = 0
      END DO
C
C-------FILL IA AND JA
      CALL PCGU7FCN(NORD,NLAY,NROW,NCOL,NIAC,NNZC,NLAPC,NUAPC,
     2                   NIAPROC,NNZAPROC,NODEC,
     3                   IAC,JAC,IUC,IXMAP,
     4                   IAPC,JAPC,IRAPC,IAPCMAP)
C
C-------SET POINTERS FOR GRID
      CALL PCGU7PSV(IGRID)
C
C-------RETURN
      RETURN
      END SUBROUTINE PCGU7AR
      
      SUBROUTINE PCGU7AP(ICNVG,KSTP,KPER,NSTP,
     &                 MXITER,KITER,
     &                 NCOL,NROW,NLAY,NODES,HDRY,IOUT,
     &                 NPC,NORD,ITER1,NITER,NNZC,NIAC,NIAPC,NIWC,
     &                 NLAPC,NUAPC,NNZAPROC,NIAPROC,
     &                 HCLOSE,RCLOSE,
     &                 HNEW,IBOUND,CR,CC,CV,HCOF,RHS,
     &                 IPRPCGU,MUTPCGU,LHCHPCGU,LRCHPCGU,IT1PCGU,
     &                 HCHGPCGU,RCHGPCGU,
     &                 NODEC,BC,XC,AC,TAPC,APC,IAC,JAC,IUC,IXMAP,IWC,
     &                 IAPC,JAPC,IRAPC,IAPCMAP,
     &                 DC,PC,QC,ZC,SCL)
C
C     ******************************************************************
C     SOLUTION BY THE CONJUGATE GRADIENT METHOD -
C                                          UP TO ITER1 ITERATIONS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(INOUT)                            :: ICNVG
      INTEGER, INTENT(IN)                               :: KSTP
      INTEGER, INTENT(IN)                               :: KPER
      INTEGER, INTENT(IN)                               :: NSTP
      INTEGER, INTENT(IN)                               :: MXITER
      INTEGER, INTENT(IN)                               :: KITER
      INTEGER, INTENT(IN)                               :: NCOL
      INTEGER, INTENT(IN)                               :: NROW
      INTEGER, INTENT(IN)                               :: NLAY
      INTEGER, INTENT(IN)                               :: NODES
      REAL, INTENT(IN)                                  :: HDRY
      INTEGER, INTENT(IN)                               :: IOUT
      INTEGER, INTENT(IN)                               :: NPC
      INTEGER, INTENT(IN)                               :: NORD
      INTEGER, INTENT(IN)                               :: ITER1
      INTEGER, INTENT(INOUT)                            :: NITER
      INTEGER, INTENT(IN)                               :: NNZC
      INTEGER, INTENT(IN)                               :: NIAC
      INTEGER, INTENT(IN)                               :: NIAPC
      INTEGER, INTENT(IN)                               :: NIWC
      INTEGER, INTENT(IN)                               :: NLAPC
      INTEGER, INTENT(IN)                               :: NUAPC
      INTEGER, INTENT(IN)                               :: NNZAPROC
      INTEGER, INTENT(IN)                               :: NIAPROC
      REAL, INTENT(IN)                                  :: HCLOSE
      REAL, INTENT(IN)                                  :: RCLOSE
      DOUBLEPRECISION, DIMENSION(NODES), INTENT(INOUT)  :: HNEW
      INTEGER, DIMENSION(NODES), INTENT(INOUT)          :: IBOUND
      REAL, DIMENSION(NODES), INTENT(IN)                :: CR
      REAL, DIMENSION(NODES), INTENT(IN)                :: CC
      REAL, DIMENSION(NODES), INTENT(IN)                :: CV
      REAL, DIMENSION(NODES), INTENT(IN)                :: HCOF
      REAL, DIMENSION(NODES), INTENT(IN)                :: RHS
      INTEGER, INTENT(IN)                               :: IPRPCGU
      INTEGER, INTENT(IN)                               :: MUTPCGU
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(INOUT) :: LHCHPCGU
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(INOUT) :: LRCHPCGU
      INTEGER, DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: IT1PCGU
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: HCHGPCGU
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: RCHGPCGU
      INTEGER, DIMENSION(NCOL,NROW,NLAY)                :: NODEC
      DOUBLEPRECISION, DIMENSION(NIAC), INTENT(INOUT)   :: BC
      DOUBLEPRECISION, DIMENSION(NIAC), INTENT(INOUT)   :: XC
      DOUBLEPRECISION, DIMENSION(NNZC), INTENT(INOUT)   :: AC
      DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: TAPC
      DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: APC
      INTEGER, DIMENSION(NIAC+1), INTENT(IN)            :: IAC
      INTEGER, DIMENSION(NNZC), INTENT(IN)              :: JAC
      INTEGER, DIMENSION(NIAC), INTENT(IN)              :: IUC
      INTEGER, DIMENSION(NIAC), INTENT(IN)              :: IXMAP
      INTEGER, DIMENSION(NIAPROC+1), INTENT(IN)         :: IAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: JAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: IRAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: IAPCMAP
C       WORKING ARRAYS
      INTEGER, DIMENSION(NIWC)          :: IWC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: R
      DOUBLEPRECISION, DIMENSION(NIAC)  :: R0
      DOUBLEPRECISION, DIMENSION(NIAC)  :: X0
      DOUBLEPRECISION, DIMENSION(NIAC)  :: DX
      DOUBLEPRECISION, DIMENSION(NIAC)  :: DC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: PC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: QC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: ZC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: SCL
C     + + + LOCAL DEFINITIONS + + +
      DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
      DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
      INTEGER :: i, j, k, n
      INTEGER :: ii
      INTEGER :: iapos
      INTEGER :: ieq
      INTEGER :: nrc
      INTEGER :: iinactive
      INTEGER :: iadiag
      INTEGER :: nrn, nrl, ncn, ncl, nln, nll
      INTEGER :: ncf, ncd, nrb, nrh, nls, nlz
      INTEGER :: ncount

      DOUBLE PRECISION :: rrhs, hhcof, rsq, fbar, fmax
      DOUBLE PRECISION :: z, b, d, e, f, h, s
      DOUBLE PRECISION :: zhnew, bhnew, dhnew, fhnew, hhnew, shnew

C     + + + FUNCTIONS + + +

C     + + + CODE + + +
      nrc     = NROW * NCOL
      iapos   = 0
      ieq     = 0
      fmax    = dzero
      ncount  = 0
      DO n = 1, NIAC
        XC(n) = HNEW(n)
        BC(n) = dzero
      END DO
      DO n = 1, NNZC
        AC(n) = dzero
      END DO
C
C-------LOOP THROUGH ALL NODES IN THE GRID AND SET UP MATRIX EQUATIONS.
C       NOTE THAT THE FORMULATION OF THESE EQUATIONS IS OPPOSITE IN SIGN
C       FROM WHAT IS GIVEN IN THE MODFLOW MANUAL SO THAT THE DIAGONAL
C       AND RHS ARE BOTH POSITIVE (LHS AND RHS ARE MULTIPLIED BY -1)
C       THIS LOOP STRUCTURE AND INDEXING IS IDENTICAL TO THAT OF PCG2 
C       AND IS BLATANTLY COPIED FROM HILL, 1990.
      LFILL: DO k = 1, NLAY
        RFILL: DO i = 1, NROW
          CFILL: DO j = 1, NCOL
C
C-------------CALCULATE 1 DIMENSIONAL SUBSCRIPT OF CURRENT CELL AND
C             INITIALIZE MATRIX COEFFICIENTS TO ZERO. CHECK IF CELL IS ACTIVE 
C             SKIP COEFFICIENT CALCULATIONS IF CELL IS INACTIVE
            n = j + (i-1) * NCOL + (k-1) * nrc
            e = dzero
            z = dzero
            b = dzero
            d = dzero
            f = dzero
            h = dzero
            s = dzero
            iinactive = 1
            ieq = ieq + 1
            IF( IBOUND(n).GT.0 ) THEN
              iinactive = 0
!              ieq       = ieq + 1
C
C---------------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR LOCATING THE 6
C               SURROUNDING CELLS
              nrn = n + NCOL
              nrl = n - NCOL
              ncn = n + 1
              ncl = n - 1
              nln = n + nrc
              nll = n - nrc
C
C---------------CALCULATE 1 DIMENSIONAL SUBSCRIPTS FOR CONDUCTANCE TO THE 6
C               SURROUNDING CELLS.
              ncf = n
              ncd = n - 1
              nrb = n - NCOL
              nrh = n
              nls = n
              nlz = n - nrc
C
C---------------STORE DOUBLE PRECISION VALUE OF RHS FOR CALCULATION OF RESIDUALS
              rrhs    =  RHS(n)
              BC(ieq) = -rrhs
C
C---------------GET CONDUCTANCES TO NEIGHBORING CELLS.  
C               ACCUMULATE CONTRIBUTIONS TO DIAGONAL COEFFICIENT. IF NEIGHBOR IS 
C               CONSTANT HEAD, MODIFY RHS AND SET OFF-DIAGONAL COEFFICIENT TO 0
              iadiag = IAC(ieq)
              iapos  = iadiag
C
C               TOP FACE
C---------------NEIGHBOR IS 1 LAYER BEHIND
              zhnew = dzero
              IF ( k.GT.1 ) THEN
                iapos = iapos + 1
                z = CV(nlz)
                e = e + z 
                zhnew = z*(HNEW(nll) - HNEW(n))
                IF ( IBOUND(nll).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -z
                END IF
                IF( IBOUND(nll).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + z*HNEW(nll) 
                  z = DZERO
                END IF
              END IF
C
C               UPPER FACE
C---------------NEIGHBOR IS 1 ROW BACK
              bhnew = dzero
              IF ( i.GT.1 ) THEN
                iapos = iapos + 1
                b = CC(nrb)
                e = e + b
                bhnew = b*(HNEW(nrl) - HNEW(n))
                IF( IBOUND(nrl).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -b
                END IF 
                IF( IBOUND(nrl).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + b*HNEW(nrl) 
                  b = dzero
                END IF 
              END IF
C
C               LEFT FACE
C---------------NEIGHBOR IS 1 COLUMN BACK
              dhnew = dzero
              IF ( j.GT.1 ) THEN
                iapos = iapos + 1
                d = CR(ncd)
                e = e + d
                dhnew = d*(HNEW(ncl) - HNEW(n))
                IF( IBOUND(ncl).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -d
                END IF
                IF( IBOUND(ncl).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + d*HNEW(ncl) 
                  d = dzero
                END IF 
              END IF
C
C               RIGHT FACE
C---------------NEIGHBOR IS 1 COLUMN AHEAD
              fhnew = dzero
              IF ( j.LT.NCOL ) THEN
                iapos = iapos + 1
                f = CR(ncf)
                e = e + f
                fhnew = f*(HNEW(ncn) - HNEW(n))
                IF( IBOUND(ncn).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -f
                END IF
                IF( IBOUND(ncn).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + f*HNEW(ncn) 
                  f = dzero 
                END IF
              END IF
C
C               LOWER FACE
C---------------NEIGHBOR IS 1 ROW AHEAD
              hhnew = dzero
              IF ( i.LT.NROW ) THEN
                iapos = iapos + 1
                h = CC(nrh)
                e = e + h
                hhnew = h*(HNEW(nrn) - HNEW(n))
                IF( IBOUND(nrn).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -h
                END IF
                IF( IBOUND(nrn).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + h*HNEW(nrn) 
                  h = dzero
                END IF 
              END IF
C
C               BOTTOM FACE
C---------------NEIGHBOR IS 1 LAYER AHEAD
              shnew = dzero
              IF ( k.LT.NLAY ) THEN
                iapos = iapos + 1
                s = CV(nls)
                e = e + s
                shnew = s*(HNEW(nln) - HNEW(n))
                IF( IBOUND(nln).GT.0 ) THEN
!                  iapos = iapos + 1
                  AC(iapos) = -s
                END IF
                IF( IBOUND(nln).LT.0 ) THEN
                  BC(ieq) = BC(ieq) + s*HNEW(nln) 
                  s = dzero
                END IF
              END IF
C
C---------------ADD HEAD-DEPENDENT BOUNDARY CONDITIONS AND STORAGE CHANGES (HCOF) TO e
              e  = e - HCOF(n)
C    
C---------------CHECK IF SURROUNDING CELLS ARE ACTIVE AND/OR HEAD-DEPENDENT BOUNDARY
C               CONDITIONS AND STORAGE CHANGES CAN OCCUR IN THE CELL (e > 0).  
C               IF SO, CALCULATE L2 NORM.  ACCUMULATE THE AVERAGE ABSOLUTE VALUE 
C               OF THE RHS VECTOR FOR ALL ACTIVE CELLS.  THIS IS USED TO SCALE THE 
C               CLOSURE CRITERIA. 
C               IF SURROUNDING CELLS ARE INACTIVE BUT CURRENT CELL IS ACTIVE, AND
C               THERE ARE NO HEAD-DEPENDENT BOUNDARY CONDITIONS AND STORAGE CHANGES
C               IN THE CURRENT CELL SET HNEW TO HDRY, IBOUND TO 0, AND CHANGE 
C               INACTIVE FLAG TO 1
              IF ( e.GT.dzero ) THEN
                hhcof  = HNEW(n)*HCOF(n)
                rsq    = rsq + (rrhs - zhnew - bhnew - dhnew - hhcof - 
     &                          fhnew - hhnew - shnew)**2
!                e      = e - HCOF(n)
                fbar   = fbar + ABS( BC(ieq) )
                ncount = ncount + 1
              ELSE
                HNEW(ieq)   = HDRY
                IBOUND(ieq) = 0
                iinactive   = 1
C                 IF INACTIVE OR CONSTANT HEAD, SET DIAGONAL TO 1.0, AND ADJUST RHS ACCORDINGLY.  
                e           = done
                BC(ieq)     = HNEW(n)
              END IF
!C    
!C---------------CHECK IF SURROUNDING CELLS ARE ACTIVE (E > 0).  IF SO, CALCULATE 
!C---------------L2 NORM.  ACCUMULATE THE AVERAGE ABSOLUTE VALUE  OF THE RHS 
!C---------------VECTOR FOR ALL ACTIVE CELLS.  THIS IS USED TO SCALE THE THE 
!C---------------CLOSURE CRITERIA. 
!C---------------IF SURROUNDING CELLS ARE INACTIVE BUT CURRENT CELL IS ACTIVE,
!C---------------SET HNEW TO HDRY, IBOUND TO 0, AND CHANGE INACTIVE FLAG TO 1
!              IF ( e.GT.dzero ) THEN
!                hhcof = HNEW(n)*HCOF(n)
!                rsq = rsq + (rrhs - zhnew - bhnew - dhnew - hhcof - 
!     &             fhnew - hhnew - shnew)**2
!                e = e - HCOF(n)
!                fbar = fbar + ABS( BC(ieq) )
!                ncount = ncount + 1
!              ELSE
!                HNEW(n)   = HDRY
!                IBOUND(n) = 0
!                iinactive = 1
!C                 IF INACTIVE OR CONSTANT HEAD, SET DIAGONAL TO 1.0, AND ADJUST RHS ACCORDINGLY.  
!                e = done
!                BC(ieq) = HNEW(n)
!              END IF
C
C-------------FIND THE MAXIMUM VALUE OF THE RHS VECTOR FOR ALL CELLS (ACTIVE
C-------------AND INACTIVE) FOR CLOSURE SCALING USED BY THE UNSTRUCTURED PCG SOLVER
              fmax = MAX( fmax, ABS( BC(ieq) ) )
C
C---------------STORE THE COEFFICENTS OF THE DIAGONAL IN A
              AC(iadiag) = e
C---------------STORE INITIAL GUESS OF HEADS
              XC(ieq) = HNEW(n)     
C
C-------------END IBOUND(N) .GT. 0
            END IF
          END DO CFILL
        END DO RFILL
      END DO LFILL
C
C-------ADD NON-ZERO VALUE TO DIAGONAL AND RIGHT-HAND SIDE IF A
C       ZERO IS PRESENT ON THE DIAGONAL
      DO n = 1, NIAC
C         SKIP INACTIVE AND CONSTANT CELLS
        IF ( IBOUND(n).LT.1 ) CYCLE
        iadiag = IAC(n)
        e      = ABS( AC(iadiag) )
        IF ( e.GT.dzero ) CYCLE
        AC(iadiag) = DONE
        BC(n)      = BC(n) + DONE
      END DO      
C
C-------SOLVE USING THE PRECONDITIONED CONJUGATE GRADIENT METHOD
      CALL SPCGUSOL(MXITER,KITER,ICNVG,NCOL,NROW,NLAY,NODES,
     &              NPC,NORD,ITER1,NITER,NNZC,NIAC,NIAPC,NIWC,
     &              NLAPC,NUAPC,NNZAPROC,NIAPROC,
     &              HCLOSE,RCLOSE,
     &              LHCHPCGU,LRCHPCGU,IT1PCGU,
     &              HCHGPCGU,RCHGPCGU,
     &              IBOUND,BC,XC,AC,TAPC,APC,IAC,JAC,IUC,IXMAP,
     &              IWC,IAPC,JAPC,IRAPC,IAPCMAP,
     &              DC,PC,QC,ZC,SCL)
C
C-------FILL HNEW WITH NEW ESTIMATE
      DO n = 1, NIAC
C         SKIP INACTIVE AND CONSTANT CELLS
        IF ( IBOUND(n).LT.1 ) CYCLE
        HNEW(IXMAP(n)) = XC(n)
!        IF ( ISNAN(XC(n)) ) THEN
!          e = XC(n)
!        ELSE
!          HNEW(IXMAP(n)) = XC(n)
!        END IF
      END DO
C
C-------IF END OF TIME STEP, PRINT # OF ITERATIONS THIS STEP
      IF ( ICNVG.NE.0 .OR. KITER.EQ.MXITER ) THEN
        IF ( MUTPCGU.LT.2 ) THEN
          IF ( KSTP.EQ.1 ) WRITE (IOUT,510)
          WRITE (IOUT,515) KITER, KSTP, KPER, NITER
C
C-----------PRINT HEAD CHANGE EACH ITERATION IF PRINTOUT INTERVAL IS REACHED
          IF ( MUTPCGU.LE.0 ) THEN
            IF ( ICNVG.EQ.0 .OR. KSTP.EQ.NSTP .OR. 
     2           MOD(KSTP,IPRPCGU).EQ.0 ) THEN
              CALL SPCG7PU(MUTPCGU,NCOL,NROW,ITER1,NITER,MXITER,IOUT,
     2                     LHCHPCGU,LRCHPCGU,IT1PCGU,HCHGPCGU,RCHGPCGU)
            END IF
          ENDIF
        ELSE IF ( MUTPCGU.EQ.3 ) THEN
          IF (ICNVG.EQ.0) THEN
              CALL SPCG7PU(0,NCOL,NROW,ITER1,NITER,MXITER,IOUT,
     2                     LHCHPCGU,LRCHPCGU,IT1PCGU,HCHGPCGU,RCHGPCGU)
          END IF
        END IF
        NITER = 0
      ENDIF
  510 FORMAT (1X,/1X)
  515 FORMAT (I6,' CALLS TO PCGU ROUTINE FOR TIME STEP',I4,
     &        ' IN STRESS PERIOD ',I4,/,I6,' TOTAL ITERATIONS')
C
C-------RETURN
      RETURN
C
      END SUBROUTINE PCGU7AP

C
C-------ROUTINE TO CALCULATE SOLVER DIMENSIONS
      SUBROUTINE PCGU7CCN(NLAY,NROW,NCOL,NIAC,NNZC,NODEC)
        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN)                              :: NLAY
        INTEGER, INTENT(IN)                              :: NROW
        INTEGER, INTENT(IN)                              :: NCOL
        INTEGER, INTENT(INOUT)                           :: NIAC
        INTEGER, INTENT(INOUT)                           :: NNZC
        INTEGER, DIMENSION(NCOL,NROW,NLAY), INTENT(INOUT):: NODEC
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k
        INTEGER :: ieq
C       + + + CODE + + +
        ieq = 0
        DO K = 1, NLAY
          DO I = 1, NROW
            DO J = 1, NCOL
              NIAC = NIAC + 1
              NNZC = NNZC + 1
              ieq  = ieq  + 1
              NODEC(J,I,K) = ieq
C               TOP FACE
              IF ( K.GT.1 ) THEN
                NNZC = NNZC + 1
              END IF
C               UPPER FACE
              IF ( I.GT.1 ) THEN
                NNZC = NNZC + 1
              END IF
C               LEFT FACE
              IF ( J.GT.1 ) THEN
                NNZC = NNZC + 1
              END IF
C               RIGHT FACE
              IF ( J.LT.NCOL ) THEN
                NNZC = NNZC + 1
              END IF
C               LOWER FACE
              IF ( I.LT.NROW ) THEN
                NNZC = NNZC + 1
              END IF
C               BOTTOM FACE
              IF ( K.LT.NLAY ) THEN
                NNZC = NNZC + 1
              END IF
            END DO
          END DO
        END DO
C
        RETURN
      END SUBROUTINE PCGU7CCN

C-------ROUTINE TO FILL IA AND JA      
      SUBROUTINE PCGU7FCN(NORD,NLAY,NROW,NCOL,NIAC,NNZC,NLAPC,NUAPC,
     2                         NIAPROC,NNZAPROC,NODEC,
     3                         IAC,JAC,IUC,IXMAP,
     4                         IAPC,JAPC,IRAPC,IAPCMAP)

        IMPLICIT NONE
C       + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN)                              :: NORD
        INTEGER, INTENT(IN)                              :: NLAY
        INTEGER, INTENT(IN)                              :: NROW
        INTEGER, INTENT(IN)                              :: NCOL
        INTEGER, INTENT(IN)                              :: NIAC
        INTEGER, INTENT(IN)                              :: NNZC
        INTEGER, INTENT(INOUT)                           :: NLAPC
        INTEGER, INTENT(INOUT)                           :: NUAPC
        INTEGER, INTENT(IN)                              :: NIAPROC
        INTEGER, INTENT(IN)                              :: NNZAPROC
        INTEGER, DIMENSION(NCOL,NROW,NLAY), INTENT(INOUT):: NODEC
        INTEGER, DIMENSION(NIAC+1), INTENT(INOUT)        :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(INOUT)          :: JAC
        INTEGER, DIMENSION(NIAC), INTENT(INOUT)          :: IUC
        INTEGER, DIMENSION(NIAC), INTENT(INOUT)          :: IXMAP
        INTEGER, DIMENSION(NIAPROC+1), INTENT(INOUT)     :: IAPC
        INTEGER, DIMENSION(NNZAPROC), INTENT(INOUT)      :: JAPC
        INTEGER, DIMENSION(NNZAPROC), INTENT(INOUT)      :: IRAPC
        INTEGER, DIMENSION(NNZAPROC), INTENT(INOUT)      :: IAPCMAP
C       + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, j, k, n
        INTEGER :: ind, ic, ieq
        INTEGER :: iz, iu, icol
        INTEGER :: nn
C       + + + CODE + + +
        ind = 0
        ic  = 0
        ieq = 0
        IALAY: DO k = 1, NLAY
          IAROW: DO i = 1, NROW
            IACOL: DO j = 1, NCOL
              ind = ind + 1
              ic  = ic + 1
              ieq = ieq + 1
              IAC(ieq) = ic
              IXMAP(ieq) = ind
              JAC(ic) = NODEC(j,i,k)
C                 TOP FACE
              if ( k.GT.1 ) then
                ic = ic + 1
                JAC(ic) = NODEC(j,i,k-1)
              END IF
C                 UPPER FACE
              IF ( i.GT.1 ) THEN
                ic = ic + 1
                JAC(ic) = NODEC(j,i-1,k)
              END IF
C                 LEFT FACE
              IF ( j.GT.1 ) THEN
                ic = ic + 1
                JAC(ic)  = NODEC(j-1,i,k)
              END IF
C                 RIGHT FACE
              IF ( j.LT.NCOL ) THEN
                ic = ic + 1
                JAC(ic) = NODEC(j+1,i,k)
              END IF
C                 LOWER FACE
              IF ( i.LT.NROW ) THEN
                ic = ic + 1
                JAC(ic) = NODEC(j,i+1,k)
              END IF
C                 BOTTOM FACE
              IF ( k.LT.NLAY ) THEN
                ic = ic + 1
                JAC(ic) = NODEC(j,i,k+1)
              END IF
            END DO IACOL
          END DO IAROW
        END DO IALAY
C
C---------SET LAST POSITION IN IAC
        IAC(NIAC+1) = ic + 1
C
C---------SET IUC
        DO n = 1, NIAC
          DO i = IAC(n), IAC(n+1)-1
            IF ( JAC(i).GT.n .AND. IUC(n).EQ.0 ) IUC(n) = i
          END DO
C           SET IUC TO FIRST ELEMENT OF THE NEXT EQUATION IF DIAGONAL
C           IS LARGEST NODE NUMBER IN AC FOR CURRENT EQUATION              
          IF ( IUC(n).EQ.0 ) IUC(n) = IAC(n+1)
        END DO
C
C---------SET INDEX FOR PRECONDITIONER IF REORDERING ILU0 AND MILU0 FOR
C         FEWER CACHE MISSES - SMITH AND ZHENG (2010)
        IF ( NORD.GT.0 ) THEN
          iz = 0
          iu = 0
C-----------LOWER
          DO n = 1, NIAC
            LOWER: DO I = IAC(n)+1, IAC(n+1)-1
              icol = JAC(i)
              IF ( icol.LT.N ) THEN
                iz = iz + 1
                JAPC(iz) = icol
                IRAPC(iz) = n
                IAPCMAP(iz) = i
              ELSE
                EXIT LOWER
              END IF
            END DO LOWER
          END DO
          NLAPC = iz
          nn    = 0
C-----------DIAGONAL AND UPPER
          DO n = NIAC, 1, -1
C             DIAGONAL
            iz = iz + 1
            iu = iu + 1
            nn = nn + 1
            IAPC(nn)    = iz
            JAPC(iz)    = n
            IRAPC(iz)   = n
            IAPCMAP(iz) = IAC(n)
C             UPPER        
            UPPER: DO i = IAC(n)+1, IAC(n+1)-1
              icol = JAC(i)
              IF ( icol.GT.n ) THEN
                iz = iz + 1
                iu = iu + 1
                JAPC(iz)    = icol
                IRAPC(iz)   = n
                IAPCMAP(iz) = i
              END IF
            END DO UPPER
          END DO
          IAPC(NIAC+1) = iz + 1
          NUAPC = iu
        END IF
C
        RETURN
      END SUBROUTINE PCGU7FCN
C
C-------ROUTINE TO SOLVE USING THE PRECONDITIONED CONJUGATE GRADIENT METHOD
      SUBROUTINE SPCGUSOL(MXITER,KITER,ICNVG,NCOL,NROW,NLAY,NODES,
     &                    NPC,NORD,ITER1,NITER,NNZC,NIAC,NIAPC,NIWC,
     &                    NLAPC,NUAPC,NNZAPROC,NIAPROC,
     &                    HCLOSE,RCLOSE,
     &                    LHCHPCGU,LRCHPCGU,IT1PCGU,
     &                    HCHGPCGU,RCHGPCGU,
     &                    IBOUND,BC,XC,AC,TAPC,APC,IAC,JAC,IUC,IXMAP,
     &                    IWC,IAPC,JAPC,IRAPC,IAPCMAP,
     &                    DC,PC,QC,ZC,SCL)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN)                               :: MXITER
      INTEGER, INTENT(IN)                               :: KITER
      INTEGER, INTENT(INOUT)                            :: ICNVG
      INTEGER, INTENT(IN)                               :: NCOL
      INTEGER, INTENT(IN)                               :: NROW
      INTEGER, INTENT(IN)                               :: NLAY
      INTEGER, INTENT(IN)                               :: NODES
      INTEGER, INTENT(IN)                               :: NPC
      INTEGER, INTENT(IN)                               :: NORD
      INTEGER, INTENT(IN)                               :: ITER1
      INTEGER, INTENT(INOUT)                            :: NITER
      INTEGER, INTENT(IN)                               :: NNZC
      INTEGER, INTENT(IN)                               :: NIAC
      INTEGER, INTENT(IN)                               :: NIAPC
      INTEGER, INTENT(IN)                               :: NIWC
      INTEGER, INTENT(IN)                               :: NLAPC
      INTEGER, INTENT(IN)                               :: NUAPC
      INTEGER, INTENT(IN)                               :: NNZAPROC
      INTEGER, INTENT(IN)                               :: NIAPROC
      REAL, INTENT(IN)                                  :: HCLOSE
      REAL, INTENT(IN)                                  :: RCLOSE
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(INOUT) :: LHCHPCGU
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(INOUT) :: LRCHPCGU
      INTEGER, DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: IT1PCGU
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: HCHGPCGU
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(INOUT) :: RCHGPCGU
      INTEGER, DIMENSION(NODES), INTENT(INOUT)          :: IBOUND
      DOUBLEPRECISION, DIMENSION(NIAC), INTENT(INOUT)   :: BC
      DOUBLEPRECISION, DIMENSION(NIAC), INTENT(INOUT)   :: XC
      DOUBLEPRECISION, DIMENSION(NNZC), INTENT(INOUT)   :: AC
      DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: TAPC
      DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: APC
      INTEGER, DIMENSION(NIAC+1), INTENT(IN)            :: IAC
      INTEGER, DIMENSION(NNZC), INTENT(IN)              :: JAC
      INTEGER, DIMENSION(NIAC), INTENT(IN)              :: IUC
      INTEGER, DIMENSION(NIAC), INTENT(IN)              :: IXMAP
      INTEGER, DIMENSION(NIAPROC+1), INTENT(IN)         :: IAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: JAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: IRAPC
      INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: IAPCMAP
C       WORKING ARRAYS
      INTEGER, DIMENSION(NIWC)          :: IWC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: R
      DOUBLEPRECISION, DIMENSION(NIAC)  :: R0
      DOUBLEPRECISION, DIMENSION(NIAC)  :: X0
      DOUBLEPRECISION, DIMENSION(NIAC)  :: DX
      DOUBLEPRECISION, DIMENSION(NIAC)  :: DC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: PC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: QC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: ZC
      DOUBLEPRECISION, DIMENSION(NIAC)  :: SCL
C     + + + LOCAL DEFINITIONS + + +
      DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
      DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
      INTEGER :: n
      INTEGER :: iiter
      INTEGER :: iicnvg
      INTEGER :: ii
      INTEGER :: ih, jh, kh, nh
      INTEGER :: ir, jr, kr, nr

      DOUBLE PRECISION :: dhclose, drclose

      DOUBLEPRECISION :: t
      DOUBLEPRECISION :: deltax, sdeltax
      DOUBLEPRECISION :: rmax, srmax
      DOUBLEPRECISION :: alpha, beta
      DOUBLEPRECISION :: rho, rho0
      
C     + + + FUNCTIONS + + +
      DOUBLEPRECISION :: SPCGUDP

C     + + + CODE + + +
      dhclose = REAL( HCLOSE, 8 )
      drclose = REAL( RCLOSE, 8 )
C
C---------UPDATE PRECONDITIONER
C         SCALE MATRIX FOR POLYNOMIAL PRECONDITIONER
        IF ( NPC.EQ.5 ) THEN
          CALL SPCGUSCL(1,NNZC,NIAC,IBOUND,AC,SCL,XC,IAC,JAC)
        END IF
        CALL SPCGUPCU(NNZC,NIAC,NIAPC,NIWC,
     2                NLAPC,NUAPC,NNZAPROC,NPC,NORD,
     3                IBOUND,AC,TAPC,APC,IAC,JAC,IUC,IWC,
     4                IAPCMAP)
C---------INITILIZE SOLUTION VARIABLE AND ARRAYS
        iiter = 0
        IF ( KITER.EQ.1 ) NITER = 0
        ICNVG  = 0
        iicnvg = 0
        alpha  = dzero
        beta   = dzero
        rho    = dzero
        rho0   = dzero
        DO n = 1, NIAC
          DC(n) = DZERO
          PC(n) = DZERO
          QC(n) = DZERO
          ZC(n) = DZERO
        END DO
C---------CALCULATE INITIAL RESIDUAL
        CALL SPCGUMV(NNZC,NIAC,IBOUND,AC,XC,DC,IAC,JAC)
        rmax = dzero
        DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE
          t     = DC(n)
          DC(n) = BC(n) - t
          t     = ABS( DC(n) )
          IF ( t.GT.rmax ) rmax  = t
        END DO
C---------CHECK IF INITIAL XC IS A SOLUTION
        IF ( rmax.EQ.dzero ) THEN
          ICNVG = 1
          GO TO 50
        END IF       
C---------INNER ITERATION          
        INNER: DO iiter = 1, ITER1
           NITER  = NITER + 1
C----------APPLY PRECONDITIONER
          SELECT CASE (NPC)
C             NO PRECONDITIONER
            CASE (1)
              DO n = 1, NIAC
                ZC(n) = DC(n)
              END DO
C             JACOBI PRECONDITIONER
            CASE (2)
              CALL SPCGUJACA(NIAC,IBOUND,APC,DC,ZC)
C             ILU0 AND MILU0 PRECONDITIONERS
            CASE (3,4)
              IF ( NORD.GT.0 ) THEN
                CALL SPCGUILU0ASZ(NNZC,NIAC,NIAPC,
     2                            NLAPC,NUAPC,NNZAPROC,NIAPROC,
     3                            IBOUND,APC,IAC,JAC,IUC,
     4                            IAPC,JAPC,IRAPC,DC,ZC)
              ELSE
                CALL SPCGUILU0A(NNZC,NIAC,NIAPC,
     2                          IBOUND,APC,IAC,JAC,IUC,DC,ZC)
              END IF
C             POLYNOMIAL PRECONDITIONER
            CASE (5)
              CALL SPCGUMV(NNZC,NIAC,IBOUND,APC,DC,ZC,IAC,JAC)
          END SELECT

          rho = SPCGUDP( NIAC, IBOUND, DC, ZC)
C-----------COMPUTE DIRECTIONAL VECTORS
          IF ( iiter.EQ.1 ) THEN
            DO n = 1, NIAC
C               SKIP INACTIVE AND CONSTANT CELLS
              IF ( IBOUND(n).LT.1 ) CYCLE
              PC(n) = ZC(n)
            END DO
          ELSE
            beta = rho / rho0
            DO n = 1, NIAC
C               SKIP INACTIVE AND CONSTANT CELLS
              IF ( IBOUND(n).LT.1 ) CYCLE
              PC(n) = ZC(n) + beta * PC(n)
            END DO
          END IF
C-----------COMPUTE ITERATES
C           UPDATE qc
          CALL SPCGUMV(NNZC,NIAC,IBOUND,AC,PC,QC,IAC,JAC)

          alpha = rho / SPCGUDP( NIAC, IBOUND, PC, QC)
C-----------UPDATE X AND RESIDUAL
          deltax = dzero
          rmax   = dzero
          DO n = 1, NIAC
C             SKIP INACTIVE AND CONSTANT CELLS
            IF ( IBOUND(n).LT.1 ) CYCLE
            t      = alpha * PC(n)
            XC(n)  = XC(n) + t
!            deltax = MAX( ABS(t), deltax )
            IF ( ABS(t).GT.deltax ) THEN
              deltax  = ABS(t)
              sdeltax = t
              nh = n
            END IF
!            IF ( ISNAN(XC(n)) ) THEN
!              t = 1.0
!            END IF
            t      = DC(n)
            t      = t - alpha * QC(n)
            DC(n)  = t
!            rmax   = MAX( ABS(t), rmax )
            IF ( ABS(t).GT.rmax ) THEN
              rmax  = ABS(t)
              srmax = t
              nr = n
            END IF
          END DO
C
C-----------STORE THE LARGEST UNSCALED HEAD CHANGE AND RESIDUAL VALUE
C           (THIS ITERATION) AND THEIR LOCATIONS.
          ii = NITER
          HCHGPCGU(ii)   = sdeltax
          CALL SPCGUCRL(NCOL,NROW,NLAY,nh,jh,ih,kh)
          LHCHPCGU(1,ii) = kh
          LHCHPCGU(2,ii) = ih
          LHCHPCGU(3,ii) = jh
C
          RCHGPCGU(ii)   = srmax
          CALL SPCGUCRL(NCOL,NROW,NLAY,nr,jr,ir,kr)
          LRCHPCGU(1,ii) = kr
          LRCHPCGU(2,ii) = ir
          LRCHPCGU(3,ii) = jr
C
          IT1PCGU(ii)    = 0
          IF ( iiter.EQ.1 ) IT1PCGU(ii) = 1
C
C           STANDARD CONVERGENCE CHECK - IDENTICAL TO PCG
          IF ( deltax.LE.dhclose .AND. rmax.LE.drclose ) iicnvg = 1
          IF ( MXITER.EQ.1 ) THEN
            IF ( iicnvg.EQ.1 ) ICNVG = 1
          ELSE
            IF ( iiter.EQ.1 .AND. iicnvg.EQ.1 ) ICNVG = 1
          ENDIF
C           IF ONLY NON-ZERO ENTRY IS ON DIAGONAL rmax WILL BE ZERO
C           THIS IS A SIMPLE SOLUTION OF ac(JA(IA(n)) x xc(n) = bc(n)
          IF ( iicnvg.LT.1 ) THEN
            IF ( rmax.EQ.DZERO ) THEN
              iicnvg = 1
              ICNVG  = 1
            END IF
          END IF
          IF ( iicnvg.EQ.1 ) EXIT INNER
C-----------SAVE CURRENT INNER ITERATES
          rho0 = rho
        END DO INNER
C
C         UNSCALE HNEW
00050   UNSCALE: IF ( NPC.EQ.5 ) THEN
          CALL SPCGUSCL(0,NNZC,NIAC,IBOUND,AC,SCL,XC,IAC,JAC)
        END IF UNSCALE
C
      RETURN
      END SUBROUTINE SPCGUSOL

      SUBROUTINE PCGU7DA(IGRID)
C  Deallocate PCGU DATA
      USE PCGUMODULE
C
      CALL PCGU7PNT(IGRID)
        DEALLOCATE(ITER1C,NPC,NORD,NITERC,NNZC,NIAC)
        DEALLOCATE(NIAPC,NIWC)
        DEALLOCATE(NLAPC,NUAPC)
        DEALLOCATE(NNZAPROC,NIAPROC)
        DEALLOCATE(HCLOSEPCGU,RCLOSEPCGU)
        DEALLOCATE(IPRPCGU,MUTPCGU)
        DEALLOCATE(LHCHPCGU)
        DEALLOCATE(LRCHPCGU)
        DEALLOCATE(IT1PCGU)
        DEALLOCATE(HCHGPCGU)
        DEALLOCATE(RCHGPCGU)
        DEALLOCATE(NODEC)
        DEALLOCATE(BC)
        DEALLOCATE(XC)
        DEALLOCATE(AC)
        DEALLOCATE(TAPC)
        DEALLOCATE(APC)
        DEALLOCATE(IAC)
        DEALLOCATE(JAC)
        DEALLOCATE(IUC)
        DEALLOCATE(IXMAP)
        DEALLOCATE(IAPC)
        DEALLOCATE(JAPC)
        DEALLOCATE(IRAPC)
        DEALLOCATE(IAPCMAP)
C       WORKING ARRAYS
        DEALLOCATE(IWC)
        DEALLOCATE(DC)
        DEALLOCATE(PC)
        DEALLOCATE(QC)
        DEALLOCATE(ZC)
        DEALLOCATE(SCL)
C
      RETURN
      END SUBROUTINE PCGU7DA
      
      SUBROUTINE PCGU7PNT(IGRID)
C  Set pointers to PCGU data for a grid
      USE PCGUMODULE
C
      ITER1C=>PCGUDAT(IGRID)%ITER1C
      NPC=>PCGUDAT(IGRID)%NPC
      NORD=>PCGUDAT(IGRID)%NORD
      NITERC=>PCGUDAT(IGRID)%NITERC
      NNZC=>PCGUDAT(IGRID)%NNZC
      NIAC=>PCGUDAT(IGRID)%NIAC
      NIAPC=>PCGUDAT(IGRID)%NIAPC
      NIWC=>PCGUDAT(IGRID)%NIWC
      NLAPC=>PCGUDAT(IGRID)%NLAPC
      NUAPC=>PCGUDAT(IGRID)%NUAPC
      NNZAPROC=>PCGUDAT(IGRID)%NNZAPROC
      NIAPROC=>PCGUDAT(IGRID)%NIAPROC
      HCLOSEPCGU=>PCGUDAT(IGRID)%HCLOSEPCGU
      RCLOSEPCGU=>PCGUDAT(IGRID)%RCLOSEPCGU
      IPRPCGU=>PCGUDAT(IGRID)%IPRPCGU
      MUTPCGU=>PCGUDAT(IGRID)%MUTPCGU
      LHCHPCGU=>PCGUDAT(IGRID)%LHCHPCGU
      LRCHPCGU=>PCGUDAT(IGRID)%LRCHPCGU
      IT1PCGU=>PCGUDAT(IGRID)%IT1PCGU
      HCHGPCGU=>PCGUDAT(IGRID)%HCHGPCGU
      RCHGPCGU=>PCGUDAT(IGRID)%RCHGPCGU
      NODEC=>PCGUDAT(IGRID)%NODEC
      BC=>PCGUDAT(IGRID)%BC
      XC=>PCGUDAT(IGRID)%XC
      AC=>PCGUDAT(IGRID)%AC
      TAPC=>PCGUDAT(IGRID)%TAPC
      APC=>PCGUDAT(IGRID)%APC
      IAC=>PCGUDAT(IGRID)%IAC
      JAC=>PCGUDAT(IGRID)%JAC
      IUC=>PCGUDAT(IGRID)%IUC
      IXMAP=>PCGUDAT(IGRID)%IXMAP
      IAPC=>PCGUDAT(IGRID)%IAPC
      JAPC=>PCGUDAT(IGRID)%JAPC
      IRAPC=>PCGUDAT(IGRID)%IRAPC
      IAPCMAP=>PCGUDAT(IGRID)%IAPCMAP
C       WORKING ARRAYS
      IWC=>PCGUDAT(IGRID)%IWC
      DC=>PCGUDAT(IGRID)%DC
      PC=>PCGUDAT(IGRID)%PC
      QC=>PCGUDAT(IGRID)%QC
      ZC=>PCGUDAT(IGRID)%ZC
      SCL=>PCGUDAT(IGRID)%SCL
C
      RETURN
      END SUBROUTINE PCGU7PNT

      SUBROUTINE PCGU7PSV(IGRID)
C  Save pointers to PCGU data
      USE PCGUMODULE
C
      PCGUDAT(IGRID)%ITER1C=>ITER1C
      PCGUDAT(IGRID)%NPC=>NPC
      PCGUDAT(IGRID)%NORD=>NORD
      PCGUDAT(IGRID)%NITERC=>NITERC
      PCGUDAT(IGRID)%NNZC=>NNZC
      PCGUDAT(IGRID)%NIAC=>NIAC
      PCGUDAT(IGRID)%NIAPC=>NIAPC
      PCGUDAT(IGRID)%NIWC=>NIWC
      PCGUDAT(IGRID)%NLAPC=>NLAPC
      PCGUDAT(IGRID)%NUAPC=>NUAPC
      PCGUDAT(IGRID)%NNZAPROC=>NNZAPROC
      PCGUDAT(IGRID)%NIAPROC=>NIAPROC
      PCGUDAT(IGRID)%HCLOSEPCGU=>HCLOSEPCGU
      PCGUDAT(IGRID)%RCLOSEPCGU=>RCLOSEPCGU
      PCGUDAT(IGRID)%IPRPCGU=>IPRPCGU
      PCGUDAT(IGRID)%MUTPCGU=>MUTPCGU
      PCGUDAT(IGRID)%LHCHPCGU=>LHCHPCGU
      PCGUDAT(IGRID)%LRCHPCGU=>LRCHPCGU
      PCGUDAT(IGRID)%IT1PCGU=>IT1PCGU
      PCGUDAT(IGRID)%HCHGPCGU=>HCHGPCGU
      PCGUDAT(IGRID)%RCHGPCGU=>RCHGPCGU
      PCGUDAT(IGRID)%NODEC=>NODEC
      PCGUDAT(IGRID)%BC=>BC
      PCGUDAT(IGRID)%XC=>XC
      PCGUDAT(IGRID)%AC=>AC
      PCGUDAT(IGRID)%TAPC=>TAPC
      PCGUDAT(IGRID)%APC=>APC
      PCGUDAT(IGRID)%IAC=>IAC
      PCGUDAT(IGRID)%JAC=>JAC
      PCGUDAT(IGRID)%IUC=>IUC
      PCGUDAT(IGRID)%IXMAP=>IXMAP
      PCGUDAT(IGRID)%IAPC=>IAPC
      PCGUDAT(IGRID)%JAPC=>JAPC
      PCGUDAT(IGRID)%IRAPC=>IRAPC
      PCGUDAT(IGRID)%IAPCMAP=>IAPCMAP
C       WORKING ARRAYS
      PCGUDAT(IGRID)%IWC=>IWC
      PCGUDAT(IGRID)%DC=>DC
      PCGUDAT(IGRID)%PC=>PC
      PCGUDAT(IGRID)%QC=>QC
      PCGUDAT(IGRID)%ZC=>ZC
      PCGUDAT(IGRID)%SCL=>SCL
C
      RETURN
      END SUBROUTINE PCGU7PSV

C
C-------SUBROUTINE TO CONVERT NODE NUMBER TO COLUMN, ROW, LAYER
      SUBROUTINE SPCGUCRL(NCOL,NROW,NLAY,INODE,J,I,K)
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN)    :: NCOL
      INTEGER, INTENT(IN)    :: NROW
      INTEGER, INTENT(IN)    :: NLAY
      INTEGER, INTENT(IN)    :: INODE
      INTEGER, INTENT(INOUT) :: J
      INTEGER, INTENT(INOUT) :: I
      INTEGER, INTENT(INOUT) :: K
C     + + + LOCAL DEFINITIONS + + +
      INTEGER :: ii, jj, kk, ij
      INTEGER :: nrc
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
      nrc = NCOL * NROW
      K   = ( INODE / nrc ) + 1
      ij  = INODE - ( K - 1 ) * nrc
      I   = ij / NCOL
      IF ( ( I * NCOL ).LT.ij ) I = I + 1
      J   = ij - ( I - 1 ) * NCOL
C---------RETURN
      RETURN
      END SUBROUTINE SPCGUCRL
      
!      SUBROUTINE SPCG7PU(HCHG,LHCH,RCHG,LRCH,ITER1,NITER,MXITER,IOUT,
!     &                   IT1,MUTPCG,NCOL,NROW)
      SUBROUTINE SPCG7PU(MUTPCGU,NCOL,NROW,ITER1,NITER,MXITER,IOUT,
     2                   LHCH,LRCH,IT1,HCHG,RCHG)
C     ******************************************************************
C     PRINT MAXIMUM HEAD CHANGE AND RESIDUAL VALUE FOR EACH ITERATION
C                           DURING A TIME STEP
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
      INTEGER, INTENT(IN) :: MUTPCGU
      INTEGER, INTENT(IN) :: NCOL
      INTEGER, INTENT(IN) :: NROW
      INTEGER, INTENT(IN) :: NITER
      INTEGER, INTENT(IN) :: MXITER
      INTEGER, INTENT(IN) :: ITER1
      INTEGER, INTENT(IN) :: IOUT
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(IN) :: LHCH
      INTEGER, DIMENSION(3,MXITER*ITER1), INTENT(IN) :: LRCH
      INTEGER, DIMENSION(MXITER*ITER1),   INTENT(IN) :: IT1
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(IN) :: HCHG
      REAL,    DIMENSION(MXITER*ITER1),   INTENT(IN) :: RCHG
C     + + + LOCAL DEFINITIONS + + +
      INTEGER :: i, j, k
      INTEGER :: jmin
      INTEGER :: l1, l2
      INTEGER :: ngrp
C     + + + FUNCTIONS + + +
C     + + + OUTPUT FORMATS + + +
2000     FORMAT(1X,/1X,'MAXIMUM HEAD CHANGE FOR EACH ITERATION',
     1       ' (1 INDICATES THE FIRST INNER ITERATION):',//
     2       1X,5('   HEAD CHANGE '),/
     3       1X,5('  LAYER,ROW,COL')/1X,75('-'))
2005          FORMAT(5(2X,I1,G12.4))
2010          FORMAT(1X,5(:'  (',I3,',',I3,',',I3,')'))
2015          FORMAT(5(4X,I1,G12.4,2X))
2020          FORMAT(1X,5(:'  (',I3,',',I5,',',I5,')'))
2025     FORMAT(1X,/1X,'MAXIMUM RESIDUAL FOR EACH ITERATION',
     1       ' (1 INDICATES THE FIRST INNER ITERATION):',//
     2       1X,5('   RESIDUAL    '),/
     3       1X,5('  LAYER,ROW,COL')/1X,75('-'))
2030     FORMAT(1X,/1X)
2035     FORMAT (1X,/1X,'MAXIMUM HEAD CHANGE FOR LAST ITER1 ITERATIONS',
     1          /1X,'(1 INDICATES THE FIRST INNER ITERATION):',//,
     2       1X,5('   HEAD CHANGE '),/
     3       1X,5('  LAYER,ROW,COL')/1X,75('-'))
2040     FORMAT(1X,/1X,'MAXIMUM RESIDUAL FOR LAST ITER1 ITERATIONS',
     1          /1X,'(1 INDICATES THE FIRST INNER ITERATION):',//,
     2       1X,5('   RESIDUAL    '),/
     3       1X,5('  LAYER,ROW,COL')/1X,75('-'))
C     + + + CODE + + +
      IF ( MUTPCGU.EQ.0 ) THEN
         WRITE(IOUT,2000)
         ngrp = (NITER-1) / 5 + 1
         DO k = 1, ngrp
            l1 = ( k - 1 ) * 5 + 1
            l2 = l1 + 4
            IF( k.EQ.ngrp ) l2 = NITER
            IF ( NCOL.LE.999 .AND. NROW.LE.999 ) THEN
              WRITE(IOUT,2005) (IT1(j),HCHG(j),j=l1,l2)
              WRITE(IOUT,2010) ((LHCH(i,j),i=1,3),j=l1,l2)
            ELSE
              WRITE(IOUT,2015) (IT1(j),HCHG(j),j=l1,l2)
              WRITE(IOUT,2020) ((LHCH(i,j),i=1,3),j=l1,l2)
            ENDIF
         END DO
         WRITE(IOUT,2025)
         DO k = 1, ngrp
            l1 = ( k - 1 ) * 5 + 1
            l2 = l1 + 4
            IF( k.EQ.ngrp ) l2 = NITER
            WRITE(IOUT,2005) (IT1(j),RCHG(j),j=l1,l2)
            WRITE(IOUT,2010) ((LRCH(i,j),i=1,3),j=l1,l2)
         END DO
         WRITE(IOUT,2030)
      ELSE
         WRITE(IOUT,2035)
         jmin = MAX( 1, NITER - ITER1 + 1 )
         ngrp = ( NITER - jmin ) / 5 + 1
         DO k = 1, ngrp
            l1 = ( k - 1 ) * 5 + jmin
            l2 = l1 + 4
            IF( k.EQ.ngrp ) l2 = NITER
            WRITE(IOUT,2005) (IT1(j),HCHG(j),j=l1,l2)
            WRITE(IOUT,2010) ((LHCH(i,j),i=1,3),j=l1,l2)
         END DO
         WRITE(IOUT,2040)
         DO k = 1, ngrp
            l1 = ( k - 1 ) * 5 + jmin
            l2 = l1 + 4
            IF( k.EQ.ngrp ) l2 = NITER
            WRITE(IOUT,2005) (IT1(j),RCHG(j),j=l1,l2)
            WRITE(IOUT,2010) ((LRCH(i,j),i=1,3),j=l1,l2)
         END DO
         WRITE(IOUT,2030)
      END IF
C---------RETURN
      RETURN
      END SUBROUTINE SPCG7PU
C
C-------ROUTINE TO SCALE THE COEFFICIENT MATRIX
      SUBROUTINE SPCGUSCL(ISCALE,NNZC,NIAC,IBOUND,AC,SCL,XC,IAC,JAC)
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: ISCALE
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(INOUT)  :: AC
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT)  :: SCL
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT)  :: XC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN)   :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)     :: JAC
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: n, i
        INTEGER :: ic
        INTEGER :: i0, i1
        DOUBLEPRECISION :: c1, c2, v
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
C
C---------SCALE THE COEFFICIENT MATRIX
        IF ( ISCALE.EQ.1 ) THEN
          DO n = 1, NIAC
C             SKIP INACTIVE AND CONSTANT CELLS
            IF ( IBOUND(n).LT.1 ) CYCLE
            ic = IAC(n)
            SCL(n) = 1.0D0 / SQRT( ABS( AC(ic) ) )
          END DO
          DO n = 1, NIAC
C             SKIP INACTIVE AND CONSTANT CELLS
            IF ( IBOUND(n).LT.1 ) CYCLE
            c1 = SCL(n)
            i0 = IAC(n)
            i1 = IAC(n+1) - 1
            DO i = i0, i1
              ic    = JAC(i)
              c2    = SCL(ic)
              v     = c1 * AC(i) * c2
              AC(i) = v
            END DO
          END DO
C---------UNSCALE THE HEAD AND THE RESIDUAL
        ELSE
          DO n = 1, NIAC
C             SKIP INACTIVE AND CONSTANT CELLS
            IF ( IBOUND(n).LT.1 ) CYCLE
            c1 = SCL(n)
            v  = XC(n)
            v = v * ( c1 * c1 )
            XC(n) = v
          END DO
        END IF
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUSCL
C
C-------ROUTINE TO UPDATE THE PRECONDITIONER
      SUBROUTINE SPCGUPCU(NNZC,NIAC,NIAPC,NIWC,
     2                    NLAPC,NUAPC,NNZAPROC,NPC,NORD,
     2                    IBOUND,AC,TAPC,APC,IAC,JAC,IUC,IWC,
     3                    IAPCMAP)
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, INTENT(IN) :: NIAPC
        INTEGER, INTENT(IN) :: NIWC
        INTEGER, INTENT(IN) :: NLAPC
        INTEGER, INTENT(IN) :: NUAPC
        INTEGER, INTENT(IN) :: NNZAPROC
        INTEGER, INTENT(IN) :: NPC
        INTEGER, INTENT(IN) :: NORD
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(IN)     :: AC
        DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: TAPC
        DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN)   :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)     :: JAC
        INTEGER, DIMENSION(NIAC), INTENT(IN)     :: IUC
        INTEGER, DIMENSION(NIWC), INTENT(IN)     :: IWC
        INTEGER, DIMENSION(NNZAPROC), INTENT(IN) :: IAPCMAP
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: n
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        SELECT CASE(NPC)
C           NO PRE-CONDITIONER
          CASE (1)
C           JACOBI PRE-CONDITIONER
          CASE (2)
            CALL SPCGUPCJ(NNZC,NIAC,IBOUND,AC,APC,IAC)
C           ILU0
          CASE (3,4)
            CALL SPCGUPCILU0(NPC,NORD,NNZC,NIAC,NIAPC,NIWC,NNZAPROC,
     2                       IBOUND,AC,TAPC,APC,IAC,JAC,IUC,IWC,IAPCMAP)
C           NEUMAN POLYNOMIAL
          CASE (5)
            CALL SPCGUNMPOL(NPC,NNZC,NIAC,NIAPC,IBOUND,AC,APC,IAC,JAC)
C           ADDITIONAL PRECONDITIONERS - ILU, etc.
        END SELECT
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUPCU
C
C-------JACOBI PRECONDITIONER - INVERSE OF DIAGONAL 
      SUBROUTINE SPCGUPCJ(NNZC,NIAC,IBOUND,AC,APC,IAC)
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(IN)      :: AC
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT)   :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN) :: IAC
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: i, n
        INTEGER :: ic0, ic1
        INTEGER :: id
        DOUBLEPRECISION :: t
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
        DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        DO n = 1, NIAC
          APC(n) = DZERO
        END DO
        DO n = 1, NIAC
C             SKIP INACTIVE AND CONSTANT CELLS
            IF ( IBOUND(n).LT.1 ) CYCLE
            id = IAC(n)
            t  = AC(id)
            IF ( ABS( t ).GT.DZERO ) t = DONE / t
            APC(n) = t
        END DO
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUPCJ

      SUBROUTINE SPCGUJACA(NIAC,IBOUND,A,D1,D2)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)    :: A
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)    :: D1
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT) :: D2
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: n
        DOUBLEPRECISION :: t, djac
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
        DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE
          t     = A(n) * D1(n)
          D2(n) = t
        END DO
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUJACA

      SUBROUTINE SPCGUPCILU0(NPC,NORD,NNZC,NIAC,NIAPC,NIWC,NNZAPROC,
     2                       IBOUND,AC,TAPC,APC,IAC,JAC,IUC,IWC,IAPCMAP)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NPC
        INTEGER, INTENT(IN) :: NORD
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, INTENT(IN) :: NIAPC
        INTEGER, INTENT(IN) :: NIWC
        INTEGER, INTENT(IN) :: NNZAPROC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(IN)     :: AC
        DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: TAPC
        DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN)   :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)     :: JAC
        INTEGER, DIMENSION(NIAC), INTENT(IN)     :: IUC
        INTEGER, DIMENSION(NIWC), INTENT(INOUT)  :: IWC
        INTEGER, DIMENSION(NNZAPROC), INTENT(IN) :: IAPCMAP
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: ic0, ic1, id0, iu1
        INTEGER :: iic0, iic1
        INTEGER :: j, n
        INTEGER :: jj, nn
        INTEGER :: jpos, jcol, jw
        INTEGER :: id
        DOUBLEPRECISION :: tl
        DOUBLEPRECISION :: t
        DOUBLEPRECISION :: rs
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
        DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        DO n = 1, NIAPC
          TAPC(n) = AC(n)
        END DO
        DO n = 1, NIAC
          IWC(n)  = 0
        END DO
        MAIN: DO n = 1, NIAC
C         SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE MAIN
          ic0 = IAC(n)
          ic1 = IAC(n+1)-1
          DO j = ic0, ic1
            jcol = JAC(j)
            IWC(jcol) = j
          END DO
          rs    = DZERO
          jpos = IAC(n)
          iu1 = IUC(n) - 1
          LOWER: DO j = ic0+1, iu1
            jpos = j
            jcol = JAC(j)
            IF ( IBOUND(jcol).LT.1 ) CYCLE LOWER
            id0 = IAC(jcol)
            tl = TAPC(j) * TAPC(id0)
            TAPC(j) = tl
            iic0 = IUC(jcol)
            iic1 = IAC(jcol+1) - 1
            DO jj = iic0, iic1
              jw = IWC(JAC(jj))
              IF ( jw.NE.0 ) THEN
                TAPC(jw) = TAPC(jw) - tl * TAPC(jj)
              ELSE
                IF ( NPC.EQ.4 ) rs = rs + tl * TAPC(jj)
              END IF
            END DO
          END DO LOWER
C           DIAGONAL - CALCULATE INVERSE OF DIAGONAL FOR SOLUTION
          id0 = IAC(n)
          tl  = TAPC(id0) - rs
          IF ( tl.GT.DZERO ) THEN
            TAPC(id0) = DONE / tl
          ELSE
            CALL USTOP('tl <= 0.0')
          END IF
C           RESET POINTER FOR IW TO ZERO
          DO j = ic0, ic1
            jcol = JAC(j)
            IWC(jcol) = 0
          END DO
        END DO MAIN
C---------MAP TAPC TO APC
        DO n = 1, NIAPC
          nn = n
          IF ( NORD.GT.0 ) nn = IAPCMAP(n)
          APC(n) = TAPC(nn)
        END DO
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUPCILU0

      SUBROUTINE SPCGUILU0A(NNZC,NIAC,NIAPC,
     2                      IBOUND,APC,IAC,JAC,IUC,R,D)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, INTENT(IN) :: NIAPC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NIAPC),  INTENT(INOUT)  :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN) :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)   :: JAC
        INTEGER, DIMENSION(NIAC), INTENT(IN)   :: IUC
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)     :: R
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT)  :: D
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: ic0, ic1, id0
        INTEGER :: jcol
        INTEGER :: j, n
        DOUBLEPRECISION :: t
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
C         FORWARD SOLVE - APC * D = R
        FORWARD: DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE FORWARD
          t   = R(n)
          ic0 = IAC(n) + 1
          ic1 = IUC(n) - 1
          LOWER: DO j = ic0, ic1
            jcol = JAC(j)
            t    = t - APC(j) * D(jcol)
          END DO LOWER
          D(n) = t
        END DO FORWARD
C         BACKWARD SOLVE - D = D / U
        BACKWARD: DO n = NIAC, 1, -1
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE BACKWARD
          id0 = IAC(n)
          ic0 = IUC(n)
          ic1 = IAC(n+1)-1
          t   = D(n)
          UPPER: DO j = ic0, ic1
            jcol = JAC(j)
            t    = t - APC(j) * D(jcol)
          END DO UPPER
C           COMPUTE D FOR DIAGONAL - D = D / U
          D(n) = APC(id0) * t
        END DO BACKWARD
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUILU0A

      SUBROUTINE SPCGUILU0ASZ(NNZC,NIAC,NIAPC,
     2                        NLAPC,NUAPC,NNZAPROC,NIAPROC,
     3                        IBOUND,
     4                        APC,IAC,JAC,IUC,
     5                        IAPC,JAPC,IRAPC,R,D)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, INTENT(IN) :: NIAPC
        INTEGER, INTENT(IN) :: NLAPC
        INTEGER, INTENT(IN) :: NUAPC
        INTEGER, INTENT(IN) :: NNZAPROC
        INTEGER, INTENT(IN) :: NIAPROC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NIAPC),  INTENT(INOUT) :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN)            :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)              :: JAC
        INTEGER, DIMENSION(NIAC), INTENT(IN)              :: IUC
        INTEGER, DIMENSION(NIAPROC+1), INTENT(IN)         :: IAPC
        INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: JAPC
        INTEGER, DIMENSION(NNZAPROC), INTENT(IN)          :: IRAPC
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)     :: R
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT)  :: D
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: ic0, ic1, id0
        INTEGER :: jcol
        INTEGER :: j, n
        INTEGER :: ic, ir, nn
        DOUBLEPRECISION :: t
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
C
C         FILL D WITH R
        FILL: DO n = 1, NIAC
          D(n) = R(n)
        END DO FILL
C         FORWARD SOLVE - APC * D = R
        FORWARD: DO n = 1, NLAPC
          ir   = IRAPC(n)
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(ir).LT.1 ) CYCLE FORWARD
          jcol = JAPC(n)
          D(ir) = D(ir) - APC(n) * D(jcol)
        END DO FORWARD
C         BACKWARD SOLVE - D = D / U
        BACKWARD: DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE BACKWARD
          id0 = IAPC(n)
          ic0 = IAPC(n) + 1
          ic1 = IAPC(n+1) - 1
          ir  = IRAPC(id0)
          t   = D(ir)
          UPPER: DO j = ic0, ic1
            jcol = JAPC(j)
            t    = t - APC(j) * D(jcol)
          END DO UPPER
C           COMPUTE D FOR DIAGONAL - D = D / U
          D(ir) = APC(id0) * t
        END DO BACKWARD
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUILU0ASZ

C-------NEUMAN POLYNOMIAL PRECONDITIONER
      SUBROUTINE SPCGUNMPOL(NPC,NNZC,NIAC,NIAPC,IBOUND,AC,APC,IAC,JAC)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NPC
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, INTENT(IN) :: NIAPC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(IN)     :: AC
        DOUBLEPRECISION, DIMENSION(NIAPC), INTENT(INOUT)  :: APC
        INTEGER, DIMENSION(NIAC+1), INTENT(IN)   :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)     :: JAC
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: ic0, ic1, id0, iu1
        INTEGER :: iic0, iic1
        INTEGER :: j, n
        INTEGER :: jj, nn
        INTEGER :: jpos, jcol, jw
        INTEGER :: id
        DOUBLEPRECISION :: tl
        DOUBLEPRECISION :: t
        DOUBLEPRECISION :: rs
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
        DOUBLEPRECISION, PARAMETER :: DONE  = 1.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUNMPOL


      SUBROUTINE SPCGUMV(NNZC,NIAC,IBOUND,A,D1,D2,IAC,JAC)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NNZC
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NNZC),  INTENT(IN)    :: A
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)    :: D1
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(INOUT) :: D2
        INTEGER, DIMENSION(NIAC+1), INTENT(IN) :: IAC
        INTEGER, DIMENSION(NNZC), INTENT(IN)   :: JAC
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: ic0, ic1
        INTEGER :: icol
        INTEGER :: m, n
        DOUBLEPRECISION :: t
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE
C               ADD DIAGONAL AND OFF-DIAGONAL TERMS
          t   = DZERO
          ic0 = IAC(n)
          ic1 = IAC(n+1)-1
          DO m = ic0, ic1
            icol = JAC(m) 
            t  = t + A(m) * D1(icol)
          END DO
          D2(n) = t
        END DO
C---------RETURN
        RETURN
      END SUBROUTINE SPCGUMV

      DOUBLEPRECISION FUNCTION SPCGUDP(NIAC,IBOUND,A,B) RESULT(C)
        IMPLICIT NONE
C     + + + DUMMY ARGUMENTS + + +
        INTEGER, INTENT(IN) :: NIAC
        INTEGER, DIMENSION(NIAC), INTENT(IN) :: IBOUND
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)    :: A
        DOUBLEPRECISION, DIMENSION(NIAC),  INTENT(IN)    :: B
C     + + + LOCAL DEFINITIONS + + +
        INTEGER :: n
        DOUBLEPRECISION, PARAMETER :: DZERO = 0.0D0
C     + + + FUNCTIONS + + +
C     + + + CODE + + +
        C = DZERO
        DO n = 1, NIAC
C           SKIP INACTIVE AND CONSTANT CELLS
          IF ( IBOUND(n).LT.1 ) CYCLE
          C = C + A(n) * B(n)
        END DO
C---------RETURN
        RETURN
      END FUNCTION SPCGUDP

